This is the in-the-weeds approach. In this document I will dig into the problem and create several functions that I can save to use later in a more streamlined document.
What do we have in the training set?
## fname label manually_verified
## 1 00044347.wav Hi-hat 0
## 2 001ca53d.wav Saxophone 1
## 3 002d256b.wav Trumpet 0
## 4 0033e230.wav Glockenspiel 1
## 5 00353774.wav Cello 1
## 6 003b91e8.wav Cello 0
The description of the data says that “a number of Freesound audio samples were automatically annotated with labels from the AudioSet Ontology … Then, a data validation process was carried out in which a number of participants did listen to the annotated sounds and manually assessed the presence/absence of an automatically assigned sound category”. And that “The non-verified annotations of the train set have a quality estimate of at least 65-70% in each category.”
With that in mind, I would prefer to use only the manually verified data.
## fname label manually_verified
## 2 001ca53d.wav Saxophone 1
## 4 0033e230.wav Glockenspiel 1
## 5 00353774.wav Cello 1
## 7 003da8e5.wav Knock 1
## 8 0048fd00.wav Gunshot_or_gunfire 1
## 11 006f2f32.wav Hi-hat 1
With that reduction, how many do we have in each category?
There are 67 sound clips labeled “Bark”.
Of course, all of these files are in the inputs folder.
Grab all of the sounds of dogs barking.
## Warning: package 'audio' was built under R version 3.4.4
## Class 'audioSample' atomic [1:620928] -3.05e-05 3.05e-05 0.00 -3.05e-05 -3.05e-05 ...
## ..- attr(*, "rate")= int 44100
## ..- attr(*, "bits")= int 16
The “rate” of 44,100 means that the pressure in front of the microphone was measured 44,100 times per second. Thus, the entire clip is 14 seconds long, and has 620,928 values.
This means we can assign each value a time in seconds.
getTimes = function(clipVals, rate=44100){
#clipVals - a vector representing a sound clip.
#rate - sample rate of the sound clip (probably 44,100)
# value => a vector representing the time measurements for those values
indecies = 1:length(clipVals)
return(indecies/rate)
}
toolkitList = c(toolkitList, "getTimes")plot(b1, type="l", col="darkblue", las=1)
start=45000
end=67000
abline(v=start, col="red"); abline(v=end, col="red")Zoom in to the range marked by the two red lines.
plot(start:end, b1[start:end], type="l", col="darkblue", las=1)
abline(v=start, col="red"); abline(v=end, col="red")A general bump in a sound wav view is referred to as an amplitude envelope. We need a way to automatically pick out these units of sound within a sound file.
## [1] 264
From how its described, that function should be just what I need, but it found an awful lot of peaks, and I don’t see an easy way to merge adjascent peaks.
So I considered another method that wasn’t necissarily made for sound but it looks like it will work.
We’ll have to flip between index as an x-axis, and time in seconds.
Pick out amplitude envelopes using the findpeaks function.
maxAmp = max(numB1)
minpeakheight = 0.1 * maxAmp # use for plotting later
threshold = 0.01 * maxAmp # use for plotting later
useFindPeaks = function(vals, sample.rate=44100,
maxAmp = max(vals),
minpeakheight = 0.1 * maxAmp,
threshold = 0.01 * maxAmp){
# vals - numeric vector, the sound values.
allEnv = findpeaks(vals, nups=5, ndowns=5,
minpeakheight = minpeakheight,
threshold = threshold,
minpeakdistance = sample.rate * .1)
allEnv = data.frame(allEnv)
names(allEnv) = c("height", "peak", "start", "end")
allEnv$peakTime = getTimes(numB1)[allEnv$peak]
return(allEnv)
}
allEnv = useFindPeaks(vals = numB1)
toolkitList = c(toolkitList, "useFindPeaks")Take a look at the results from findpeaks.
xlim=NULL
plot(timesB1, numB1, type="l", col="darkblue", las=1, xlim=xlim)
points(allEnv$peakTime, allEnv[,"height"], col="red", pch=16)
abline(v=getTimes(numB1)[allEnv$start], col="green")
abline(v=getTimes(numB1)[allEnv$end], col="red")xlim=c(1.10, 1.5)
plot(timesB1, numB1, type="l", col="darkblue", las=1, xlim=xlim)
points(allEnv$peakTime, allEnv[,"height"], col="red", pch=16)
abline(v=getTimes(numB1)[allEnv$start], col="green")
abline(v=getTimes(numB1)[allEnv$end], col="red")I was able to get reasonable peaks with findpeaks, but I didn’t like the start and end that was assigned to each peak. I want to use the peaks from findpeaks, but define the start and stop of each peak based on when the sound returns to some threshold.
To do this, I want to look for continuous areas of low values.
A value is only true if the surrounding n values are also true.
noLonelyTrue = function(bool, n=1, left=n, right=n, naEdges=NA){
# bool - a vectore of boolean values
# n - how many values on either side must be true for an original value to still be true
# left, right - number of values to the left (preceeding) or right (after) that must be true
end = length(bool)
# matrices are faster than loops.
# any given row in this matrix can be read as "is the original value true, how about the one before it, and before that..."
# the offset by one is acheived by adding a value at the end.
toTheLeft = matrix(data=rep(c(bool,NA),left+1)[1:(end * (left+1))], nrow=length(bool), ncol=left+1, byrow = F)
# same as above, but the other direction
toTheRight = matrix(data=rep(c(rev(bool), NA), right+1)[1:(end * (right+1))], nrow=length(bool), ncol=right+1, byrow = F)
# true values are treated as 1 in rowSums
newBool = rowSums(toTheLeft)==left+1 & rev(rowSums(toTheRight))==right+1
# the first and last n values don't have enough neighbors, by default they are NA values.
if (!is.na(naEdges)){
newBool[1:left] = naEdges
newBool[(length(newBool)-right+1):length(newBool)] = naEdges
}
return(newBool)
}
bool = c(T,T,F,T,T,T,T,T,T,T,T,F,F,T,T,T,T)
noLonelyTrue(bool, n=1)## [1] NA FALSE FALSE FALSE TRUE TRUE TRUE TRUE TRUE TRUE FALSE
## [12] FALSE FALSE FALSE TRUE TRUE NA
Now define start and stop of peask based on when there is a continuous low value.
plot(timesB1, numB1, type="l", col="darkblue", las=1, xlim=xlim)
points(allEnv$peakTime, allEnv[,"height"], col="red", pch=16)
abline(h=threshold, col="pink")
getPeaksEdges <- function(vals, peaks, threshold, minNumBelow=20){
# vals - numeric vector, the sound values.
# peaks - the indecies for the peaks
# threshold - scaler, below this value is considered quite, a pause.
# minNumBelow - integer, number of consecutive values below threshold to be called a pause.
######
# Make a matrix where each position in the sound is a row and each peak is a column.
# Make a boolean matrix describing weather each position is below threshold.
# the signal may dip below the threshold, only consider it a pause if it stays below threshold
nrow=length(vals)
ncol=length(peaks)
bt1 = vals < threshold
bt2 = noLonelyTrue(bt1, n=minNumBelow, naEdges = F)
bt = matrix(data=bt2, byrow = F, nrow=nrow, ncol=ncol)
# Make a matrix with the same layout describing weather the value (row) is before a give peak (column).
beforePeak = rep(1:length(vals), length(peaks)) < rep(peaks, each=length(vals))
bp = matrix(data=beforePeak, byrow=F, nrow=nrow, ncol=ncol)
# Combine these to get values that are pause areas before the peak.
# For each column (peak) keep the true value with highest posible index, that is the peak start.
isBefore = bt & bp
startInd = apply(isBefore, 2, function(x){max(which(x))})
# Use the mirror image of the process above to get the end for each peak.
isAfter = bt & !bp
endInd = apply(isAfter, 2, function(x){min(which(x))})
return(data.frame(peak=peaks, startInd=startInd, endInd=endInd))
}
toolkitList = c(toolkitList, "getPeaksEdges")
peakEnds = getPeaksEdges(numB1, peaks=allEnv$peak, threshold=threshold)
abline(v=getTimes(numB1)[peakEnds$startInd], col="green")
abline(v=getTimes(numB1)[peakEnds$endInd], col="red")With this method, there are some peaks that are close enough to each other that they have the same start and end postions. I don’t want to copy the same slice and mistake it for a repeating pattern, so I’m going to remove those duplicates.
allEnv$start = peakEnds$startInd
allEnv$end = peakEnds$endInd
range = paste(allEnv$start, allEnv$endInd)
allEnv = allEnv[!duplicated(range),]xlim=NULL
plot(timesB1, numB1, type="l", col="darkblue", las=1, xlim=xlim)
points(timesB1[allEnv$peak], allEnv[,"height"], col="red", pch=16)
abline(h=threshold, col="pink")
abline(v=timesB1[allEnv$start], col="green")
abline(v=timesB1[allEnv$end], col="red")This gives us 20 amplitude envelopes in this file, and most of them look reasonable.
Lets pretend we sample some really ideal sound at a rate of 10 measurments per second. Let x be your time line in seconds.
example.rate = 10
x <- seq(-30, 70, by = 1/example.rate)
wav.freq = c(0.8, 1.5, 2)
wav.1 = sin(wav.freq[1] * 2 * pi * x)
wav.2 = sin(wav.freq[2] * 2 * pi * x)
wav.3 = rep(0, length(x))
xmid = seq(0, 20, by = 1/example.rate)
someMidXs = which(x>0 & x < max(x)/3)
xmid = x[someMidXs]
wav.3[someMidXs] = 1.6 * sin(wav.freq[3] * 2 * pi * xmid)
wav.sum = wav.1 + wav.2 + wav.3
waves = list("wave1" = wav.1,
"wave2" = wav.2,
"wave3" = wav.3,
"sum" = wav.sum)The Forier transfrom is beautifully explained by 3blue1brown in his video “But what is the Fourier Transform? A visual introduction.” I’m using that as my bases for making the following function.
fourtrans = function(winding, g, t){
# winding - scaler, the "winding frequency"
# g - a vector describing the original wave.
# t - a vector of times corresponding to g; g and t must be the same length
i = complex(imaginary=1)
mass = g * exp( -2 * pi * i * winding * t)
return(sum(sum(mass))) # the double sum is diliberate.
}
toolkitList = c(toolkitList, "fourtrans")Do the fourier transform on the example waves.
ex.freqs = seq(0.001, 3, 0.001)
sum.ft = sapply(ex.freqs, fourtrans, g=wav.sum, t=x)
plot(ex.freqs, abs(sum.ft), type="l", las=1, xlab="frequency")
polygon(x=c(max(ex.freqs), min(ex.freqs), ex.freqs), c(0,0, abs(sum.ft)), col="gray")
for (i in 1:length(wav.freq)){
abline(v=wav.freq[i], col=wav.colors[i])
text(x=wav.freq[i], y=par("usr")[4], xpd=T, pos=3, col=wav.colors[i], labels=wav.freq[i])
}I see what I want. I see peaks that corresond to my original frequencies. But how to extract them when I don’t know the values and I don’t know how many there are? If I take every peak, I will have too many.
getMaxima = function(vector){
n=length(vector)
left = vector[1:(n-2)]
main = vector[2:(n-1)]
right = vector [3:n]
indecies = which(main > left & main > right)
# these indeces are for 'main' which is shifted from the original by 1.
return(indecies + 1)
}
rawMaxima = getMaxima(abs(sum.ft))
toolkitList = c(toolkitList, "getMaxima")Using the raw values and taking EVERY maximum gives me 227 frequencies.
plot(ex.freqs, abs(sum.ft), type="l", las=1, xlab="frequency")
polygon(x=c(max(ex.freqs), min(ex.freqs), ex.freqs), c(0,0, abs(sum.ft)), col="gray")
points(col="red", pch=16, x=ex.freqs[rawMaxima], abs(sum.ft)[rawMaxima])
legend(legend=length(rawMaxima), x="topright", text.col="red", pch=16, col="red")I could order them by their intensity (y-axis), but without knowing how many there are, I might stop with only 2 or take on several too many.
I know I can count on having a lot of little hills that don’t matter, so I’m going take the 75th percentile and use that a bench mark for “surely small”. Any maxima that are no higher than that, can be considered noise.
dotColor = rep("red", length(rawMaxima))
dotColor[abs(sum.ft)[rawMaxima] <= small] = "lightblue"
plot(ex.freqs, abs(sum.ft), type="l", las=1, xlab="frequency")
polygon(x=c(max(ex.freqs), min(ex.freqs), ex.freqs), c(0,0, abs(sum.ft)), col="gray")
points(col=dotColor, pch=16, x=ex.freqs[rawMaxima], abs(sum.ft)[rawMaxima])
abline(h=small, lty=2, col="lightblue")
legend(legend=sum(dotColor=="red"), x="topright", text.col="red", pch=16, col="red")Smooth out the values to remove unimportant peaks.
n.panels = 8
par(mfrow=c(n.panels,1), las=1, mar=c(0,3,0,0), oma=c(4,0,3,0))
plotSmoothMaxes = function(x, y, dotsAt, thresh){
#x - the x values of points to plot
#y - the y values of points to plot
#dotsAt - the indices of points to highlight (the maxes)
#thresh - the threshold y-value below which to ignore the maxes
plot(x, y, type="l", las=1, xlab="frequency", xaxt="n")
polygon(x=c(max(x), min(x), x), c(0,0, y), col="gray")
dotCol = rep("red", length(dotsAt))
dotCol[y[dotsAt] <= thresh] = "lightblue"
abline(h=thresh, lty=2, col="lightblue")
points(col=dotCol, pch=16, x=x[dotsAt], y[dotsAt])
}
toolkitList = c(toolkitList, "plotSmoothMaxes")
span = seq(0.01, .2, length.out = 40)
n.span = length(span)
drawThese = round(seq(1, n.span, length.out = n.panels))
drawList = list()
numMaxes = rep(NA, length(span))
for (i in 1:n.span){
ft.lo = loess(abs(sum.ft) ~ ex.freqs, span=span[i])
smoother = predict(ft.lo, ex.freqs)
maxes = getMaxima(smoother)
numMaxes[i] = sum(smoother[maxes] > small)
if (i %in% drawThese) {
plotSmoothMaxes(x=ex.freqs, y=smoother, dotsAt = maxes, thresh=small)
abline(v=wav.freq, col=wav.colors)
legend(x="topright", legend=paste("span =", span[i]), bty="n")
}
}
title("smoothed data with different 'span'", outer = T)Using a span of 0.07 is just enough to only get the main peaks. We still get the right number of peaks even with a span of more than twice that. As we increase the span, we see the values we would get from our maxima shift away from the accurate values (the red dots move away from the red, green and gold lines).
Now we need a rule a computer can follow that will help it come to similar conclusions about optimizing the span for a given clip. There is probably some really good algorithm for doing that, but for now this will do:
# What is the most common number of peaks (across different span parameters)
# That must be right. What is the lowest span that gives me that many? That must be best.
getBestSpan = function(numMaxes, span, showplot=F){
tb = table(numMaxes)
tb = tb[order(tb, decreasing = T)]
mode = as.numeric(names(tb)[1])
w = which(numMaxes == mode)
lowSpan = span[min(w)]
if (showplot){
plot(span, numMaxes, type = "l", las=1)
abline(v=lowSpan, col="red")
abline(h=mode, col="blue", lty=2)
legend(x="topright", legend=paste("span =", lowSpan), bty="n")
title("Optimal Span")
}
return(lowSpan)
}
toolkitList = c(toolkitList, "getBestSpan")
lowSpan = getBestSpan(numMaxes, span)Based on this logic, 0.068 is the optimum span parameter for this clip.
getBestFreqs = function(ft, freqs, lowSpan, small, showplot=F){
ft.lo.best = loess(abs(ft) ~ freqs, span=lowSpan)
smoother = predict(ft.lo.best, freqs)
maxes = getMaxima(smoother)
maxes = maxes[smoother[maxes] > small]
freq.of.interest = freqs[maxes]
if(showplot){
plotSmoothMaxes(x=freqs, y=smoother, dotsAt = maxes, thresh=small)
}
return(freq.of.interest)
}
freq.of.interest = getBestFreqs(sum.ft, ex.freqs, lowSpan, small=small)And based on that value for span, the frequencies of interest are:
## [1] 0.798 1.500 1.997
The peak for the wave 3 frequency is shorter and wider than the other two. This makes sense since it was only present in a small portion of the clip while the other two were present consitently. Lets apply a sliding window approach.
Define the windows.
winSize = 100 # how may values wide each window is
shift = 10 # how many values from the start of one window to the start of the next
len = length(wav.sum) # how many total values in the clip
winStarts = seq(1, len-winSize, shift)
winEnds = seq(winSize, len, shift)
nw = length(winStarts) # number of windowsFor each window, do the fourier transform and save the value for each frequency.
#ex.freqs = seq(0.001, 1, 0.001)
shiftingFT = matrix(NA, ncol = length(ex.freqs), nrow=nw,
dimnames=list(startTime=x[winStarts], frequency=ex.freqs))
for (i in 1:nw){
win = wav.sum[winStarts[i]:winEnds[i]]
shiftingFT[i,] = sapply(ex.freqs, fourtrans, g=win, t=x[winStarts[i]:winEnds[i]])
}
shiftingFT.real = abs(shiftingFT)midTime=round((x[winStarts]+x[winEnds])/2)
image(x=midTime, y=ex.freqs, z=shiftingFT.real, las=1, col=heat.colors(length(shiftingFT.real)))For just the frequences of intersest, plot their intensity over time (across the sliding windows).
freq.of.interest = freq.of.interest [order(freq.of.interest)] # make sure they are in the original order
freqName = as.character(freq.of.interest)
plot(1,1, xlim=range(x), ylim=c(0,300), type="n",
xlab="time", ylab="frequency intensity", las=1)
for (i in 1:3){
fq = freqName[i]
lines(x=midTime, y=shiftingFT.real[,fq], col=wav.colors[i])
}
legend(x="topright", legend = round(freq.of.interest,3), col=wav.colors, lty=1)The only thing we did to optimize the window size and shift size was… well nothing. I looked at the picture and adjusted it. I don’t know any good rules there, except hope that the same parameters work well over most cases.
Given this information, can you get the original data back? Because my measurements are for time intervales, and each interval must be >0, I don’t have values for the begining and end. We loose the edges of the data.
rebuilt = matrix(ncol=length(freq.of.interest),
nrow=length(x),
dimnames=list(x=x, freq=freqName))
for (i in 1:length(freq.of.interest)){
fq = freqName[i]
#.1 was the lowest span that I tried that didn't produce an error
intensityPerWindow.lo = loess(shiftingFT.real[,fq] ~ midTime, span=0.1)
amplitudeFactor = predict(intensityPerWindow.lo, x)
rebuilt[,i] = sin(freq.of.interest[i] * 2 * pi * x) * amplitudeFactor
}Normalize the amplitude to roughly match the original.
Not perfect… but not too bad. I think real sound data will have more cycles per unit, and I think that will mean an increase in resolution. It would also look cleaner to have used a cutoff rather than scaling the amplitude by the intensity, but I don’t think that is as true for real data as it is for this simulation.
Here’s the original for comparison:
The original wave (g) is represented by our input data. For t, we need to create a vector of time measurements corresponding to the values in our sound clip.
Based on Wikipedia’s Audio Frequency page, sound frequencies that we can hear will be covered by a range of about 20 to 20,000 Hz.
Do the Fourier Transform for the arf sample.
# plot(freq, abs(arf.ft), type="l", las=1, xlab="frequency")
# title("Fourier transformed dog bark")
# polygon(x=c(max(freq), min(freq), freq), c(0,0, abs(arf.ft)), col="gray")
# points(col="red", pch=16, x=freq[arfMaxes], y=abs(arf.ft)[arfMaxes])
plotSmoothMaxes(x=freq, y=abs(arf.ft),
dotsAt = arfMaxes,
thresh=quantile(abs(arf.ft)[getMaxima(abs(arf.ft))], .95))Which frequencies are of interest?
# Earlier I was interested in a lot of plotting.
# Here, just put the important bits into a function, no plots
getFreqOfInterest.old = function(ft, freqs,
small = quantile(abs(ft)[getMaxima(abs(ft))], .95),
span = seq(0.01, .2, length.out = 40)){
# ft = fourier transformed data
# ex.freqs - frequences to iterate over
# small = frequencies with intensities lower than this are discarded
# span - a vector of spans to iterate through
n.span = length(span)
numMaxes = rep(NA, length(span))
for (i in 1:n.span){
tryCatch(expr = {
numMaxes[i] = suppressWarnings(expr={
length(getBestFreqs(ft, freqs, lowSpan=span[i], small=small, showplot = T))
})
}, finally = {
#rm(ft.lo); rm(smoother); rm(maxes)
})
}
lowSpan = getBestSpan(numMaxes, span, showplot = T)
freq.of.interest = getBestFreqs(ft, freqs, lowSpan, small=small, showplot = T)
return(freq.of.interest)
}
getPeakFreqs <- function(ft, small=NULL){
if (is.null(small)){
fracMax = max(ft)/10
quant = quantile(abs(ft)[getMaxima(abs(ft))], .95)
small = max(fracMax, quant)
}
return(findpeaks(abs(ft), minpeakheight = small,
nups = 2, ndowns = 2, minpeakdistance = 5))
}
freqPeaks = getPeakFreqs(ft=abs(arf.ft))
arfFOI = freq[freqPeaks[,2]]
toolkitList = c(toolkitList, "getPeakFreqs", "getFreqOfInterest.old")
#arfFOI = getFreqOfInterest(ft=arf.ft, freqs=freq)## [1] 860 500
For each window, do the fourier transform and save the value for each frequency.
getShiftingFT = function(rawSound, sampleRate,
winSize = 5000, shift = 1000,
freqs=freq){
# rawSound - values from a sound file
# sampleRate - number of values per second in rawSound
# winSize - number of values per window
# shift - number of values to shift over to start the next window
# freq - freqencies to measure
len = length(rawSound) # how many total values in the clip
times = getTimes(rawSound, rate = sampleRate)
winStarts = seq(1, len-winSize, shift)
winEnds = seq(winSize, len, shift)
midTime = (times[winStarts] + times[winEnds])/2
nw = length(winStarts) # number of windows
shiftingFT = matrix(NA, ncol = length(freqs), nrow=nw,
dimnames=list(midTime=round(midTime,3), frequency=freqs))
for (i in 1:nw){
win = rawSound[winStarts[i]:winEnds[i]]
shiftingFT[i,] = sapply(freqs, fourtrans, g=win, t=times[winStarts[i]:winEnds[i]])
}
shiftingFT.real = abs(shiftingFT)
return(shiftingFT.real)
}
toolkitList = c(toolkitList, "getShiftingFT")
arfMatrix = getShiftingFT(rawSound = arf, sampleRate = sample.rate)midTime = as.numeric(rownames(arfMatrix))
image(x=midTime, y=freq, z=arfMatrix, las=1, col=heat.colors(length(arfMatrix)))Plotting all of the frequencies shows a lot of blank space, just show the freqencies up to 2x the highest freq of interest.
maxFreqShow = max(arfFOI) * 1.2
w = which(freq <= maxFreqShow)
plotMatrix = arfMatrix[,w]
image(x=midTime, y=as.numeric(colnames(plotMatrix)), z=plotMatrix, las=1, col=heat.colors(length(plotMatrix)), ylab="frequency (Hz)")#freq.of.interest = freq.of.interest [order(freq.of.interest)] # make sure they are in the original order
freqName = as.character(arfFOI)
arf.col = wav.colors # they can just happen to be the same
plot(1,1, xlim=range(midTime), ylim=c(0,max(arfMatrix)), type="n",
xlab="time", ylab="frequency intensity", las=1)
for (i in 1:length(arfFOI)){
fn = freqName[i]
lines(x=midTime, y=arfMatrix[,fn], col=arf.col[i])
}
legend(x="topright", legend = round(arfFOI,3), col=arf.col, lty=1)Hear each frequency alone.
for (i in 1:length(arfFOI)){
toneTime = c(times, times+max(times), times+(2*max(times)))
s = sin(arfFOI[i] * 2 * pi * toneTime)
tone = audioSample(s, rate=sample.rate, bits=16)
play(tone)
wait(2)
}Given this information, can you get the original data back? Because my measurements are for time intervales, and each interval must be >0, I don’t have values for the begining and end. We loose the edges of the data.
rebuilt = matrix(ncol=length(arfFOI),
nrow=length(times),
dimnames=list(time=times, freq=freqName))
for (i in 1:length(arfFOI)){
fq = freqName[i]
intensityPerWindow.lo = loess(c(0,arfMatrix[,fq],0) ~ c(min(times),midTime,max(times)))
amplitudeFactor = predict(intensityPerWindow.lo, times)
rebuilt[,i] = sin(arfFOI[i] * 2 * pi * times) * amplitudeFactor
}Normalize the amplitude to roughly match the original.
Smothing by maxes–this is something we might explore as a way to summarize an overal shape.
## R version 3.4.1 (2017-06-30)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS Sierra 10.12.6
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] audio_0.1-5.1 seewave_2.1.0 pracma_2.1.4
##
## loaded via a namespace (and not attached):
## [1] Rcpp_0.12.17 digest_0.6.15 rprojroot_1.3-2 MASS_7.3-50
## [5] backports_1.1.2 signal_0.7-6 magrittr_1.5 evaluate_0.10.1
## [9] stringi_1.2.3 rmarkdown_1.10 tools_3.4.1 tuneR_1.3.3
## [13] stringr_1.3.1 yaml_2.1.19 compiler_3.4.1 htmltools_0.3.6
## [17] knitr_1.20